#include<iostream>
#include<cmath>
#include <vector>
#include <fstream>


// 递归计算差商的函数，会爆内存
double get_diff_quotient(const std::vector<double> f_values, const std::vector<double> x_values, int n) {
    if (n == 0){return f_values[0];}
    else
    {
        std::vector<double> x_values_pop = x_values;
        std::vector<double> f_values_pop = f_values;
        x_values_pop.erase(x_values_pop.begin());
        f_values_pop.erase(f_values_pop.begin());
        return (get_diff_quotient(f_values_pop, x_values_pop, n-1) - get_diff_quotient(f_values, x_values, n-1)) / (x_values[n] - x_values[0]);
    }
}
// 动态规划计算差商的改进函数
void compute_diff_quotients(const std::vector<double>& x_values, const std::vector<double>& f_values, std::vector<std::vector<double>>& diff_quotients) {
    int n = x_values.size();
    // 填充差商表的第一列
    for (int i = 0; i < n; i++) {
        diff_quotients[i][0] = f_values[i];
    }

    // 使用动态规划方法填充差商表的其余部分
    for (int j = 1; j < n; j++) {
        for (int i = j; i < n; i++) {
            diff_quotients[i][j] = (diff_quotients[i][j-1] - diff_quotients[i-1][j-1]) / (x_values[i] - x_values[i-j]);
        }
    }
}
/*
double diff_quotient(double f[], double x[], int n) {
    if (n == 0)
        return f[0];
    else
        return (diff_quotient(f+1, x+1, n-1) - diff_quotient(f, x, n-1)) / (x[n-1] - x[0]);
}
*/

void multiplyPolynomials(const std::vector<double>& poly1, const std::vector<double>& poly2, std::vector<double>& polyout) {
    int size1 = poly1.size();
    int size2 = poly2.size();
    //std::cout << "multply.size1: " << size1 << "multply.size2: " << size2 << std::endl;
    int sizeout = (size1 + size2 - 1);
    // 初始化结果多项式的系数向量
    polyout.resize(sizeout,0);
    // 计算结果多项式的每一项系数
    for (int i = 0; i < size1; i++) {
        for (int j = 0; j < size2; j++) {
            polyout[i + j] += poly1[i] * poly2[j];
        }
    }
    polyout.resize(sizeout);
    return ;
}
void addPolynomials(const std::vector<double>& poly1, const std::vector<double>& poly2, std::vector<double>& polyout) {
    int size1 = poly1.size();
    int size2 = poly2.size();
    //std::cout << "add.size1: " << size1 << "add.size2: " << size2 << std::endl;
    int sizeout = std::max(size1, size2);
    // 初始化结果多项式的系数向量
    polyout.resize(sizeout,0);
    // 计算结果多项式的每一项系数
    for (int i = 0; i < size1; i++) {
            polyout[i] += poly1[i];
    }
    for (int i = 0; i < size2; i++) {
            polyout[i] += poly2[i];
    }
    polyout.resize(sizeout);
    return;
}

// 读取数据文件的函数
void readDataFile(const std::string& filename, std::vector<double>& x_values, std::vector<double>& f_values) {
    // 打开文件
    std::ifstream file(filename);
    
    if (!file) {
        std::cout << "无法打开文件 " << filename << std::endl;
        return;
    }

    // 读取数据
    double x, f;
    while (file >> x >> f) {
        x_values.push_back(x);
        f_values.push_back(f);
    }

    // 关闭文件
    file.close();
}

int main() {
    // 文件名
    std::string filename = "data_no3.txt";

    // x值和f(x)值的向量
    std::vector<double> x_values;
    std::vector<double> f_values;
    std::vector<double> output_para(10,0);
    std::vector<double> output_para_temp(1,0);
    std::vector<double> output_para_temp2(1,1);
    std::vector<double> multply_poly(2,0);

    int n;
    int i,j;
    //double diff_quotient;
    // 调用函数读取数据
    readDataFile(filename, x_values, f_values);
    n = x_values.size();
    // 初始化差商表
    std::vector<std::vector<double>> diff_quotients(n, std::vector<double>(n, 0));
    compute_diff_quotients(x_values,f_values,diff_quotients);
    // 输出读取的数据
    for (i = 0; i < n; i++) {
        std::cout << "x: " << x_values[i] << ", f(x): " << f_values[i] <<"i = "<< i << std::endl;
    }
    std::cout <<"n = "<< n << std::endl;
    for(i=0;i<n;i++){
        std::cout<<"get diff"<< std::endl;
        //diff_quotient=get_diff_quotient(x_values,f_values,i);
        std::cout <<"diff_quotient = "<< diff_quotients[i][i];
        std::cout<<"multipy"<< std::endl;
        if(i!=0){
            multply_poly[0] = 1;
            multply_poly[1] = -1*f_values[i-1];
            multiplyPolynomials(output_para_temp2,multply_poly,output_para_temp);
            output_para_temp2 = output_para_temp;
            }
        std::cout<<"add"<< std::endl;
        for(j=0;j<=i;j++){
        output_para[j] += diff_quotients[i][i]* output_para_temp[j];
        }
    }
    for (i = 0; i < n; i++) {
        std::cout << i << "次项系数: " << output_para[i] << std::endl;
    }
    std::ofstream outfile("out.txt");
    for (i = 0; i < n; i++) {
        outfile << i << "次项系数: " << output_para[i] << std::endl;
    }
    
    outfile.close();

    return 0;
}
